home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PIBSIGS / ALGAMA.PAS next >
Pascal/Delphi Source File  |  1986-06-22  |  9KB  |  248 lines
  1. (*-------------------------------------------------------------------------*)
  2. (*               ALGama  -- Logarithm of Gamma Distribution                *)
  3. (*-------------------------------------------------------------------------*)
  4.  
  5. FUNCTION ALGama( Arg : REAL ) : REAL;
  6.  
  7. (*-------------------------------------------------------------------------*)
  8. (*                                                                         *)
  9. (*       Function:  ALGama                                                 *)
  10. (*                                                                         *)
  11. (*       Purpose:  Calculates Log (base E) of Gamma function               *)
  12. (*                                                                         *)
  13. (*       Calling Sequence:                                                 *)
  14. (*                                                                         *)
  15. (*            Val := ALGama( Arg )                                         *)
  16. (*                                                                         *)
  17. (*                 Arg   --- Gamma distribution parameter (Input)          *)
  18. (*                 Val   --- output Log Gamma value                        *)
  19. (*                                                                         *)
  20. (*       Calls:  None                                                      *)
  21. (*                                                                         *)
  22. (*       Called By:  Many (CDBeta, etc.)                                   *)
  23. (*                                                                         *)
  24. (*       Remarks:                                                          *)
  25. (*                                                                         *)
  26. (*            Minimax polynomial approximations are used over the          *)
  27. (*            intervals [-inf,0], [0,.5], [.5,1.5], [1.5,4.0],             *)
  28. (*            [4.0,12.0], [12.0,+inf].                                     *)
  29. (*                                                                         *)
  30. (*            See Hart et al, "Computer Approximations",                   *)
  31. (*            Wiley(1968), p. 130F, and also                               *)
  32. (*                                                                         *)
  33. (*            Cody and Hillstrom, "Chebyshev approximations for            *)
  34. (*            the natural logarithm of the Gamma function",                *)
  35. (*            Mathematics of Computation, 21, April, 1967, P. 198F.        *)
  36. (*                                                                         *)
  37. (*                                                                         *)
  38. (*            There are some machine-dependent constants used --           *)
  39. (*                                                                         *)
  40. (*               Rmax   --- Largest value for which ALGama                 *)
  41. (*                          can be safely computed.                        *)
  42. (*               Rinf   --- Largest floating-point number.                 *)
  43. (*               Zeta   --- Smallest floating-point number                 *)
  44. (*                          such that (1 + Zeta) = 1.                      *)
  45. (*                                                                         *)
  46. (*-------------------------------------------------------------------------*)
  47.  
  48. VAR
  49.    Rarg        : REAL;
  50.    Alinc       : REAL;
  51.    Scale       : REAL;
  52.    Top         : REAL;
  53.    Bot         : REAL;
  54.    Frac        : REAL;
  55.    Algval      : REAL;
  56.  
  57.    I           : INTEGER;
  58.    Iapprox     : INTEGER;
  59.    Iof         : INTEGER;
  60.    Ilo         : INTEGER;
  61.    Ihi         : INTEGER;
  62.  
  63.    Qminus      : BOOLEAN;
  64.    Qdoit       : BOOLEAN;
  65.  
  66.  
  67. (* Structured *) CONST
  68.  
  69.    P           : ARRAY [ 1 .. 29 ] OF REAL =
  70.                  ( 4.12084318584770E+00 ,
  71.                    8.56898206283132E+01 ,
  72.                    2.43175243524421E+02 ,
  73.                   -2.61721858385614E+02 ,
  74.                   -9.22261372880152E+02 ,
  75.                   -5.17638349802321E+02 ,
  76.                   -7.74106407133295E+01 ,
  77.                   -2.20884399721618E+00 ,
  78.  
  79.                    5.15505761764082E+00 ,
  80.                    3.77510679797217E+02 ,
  81.                    5.26898325591498E+03 ,
  82.                    1.95536055406304E+04 ,
  83.                    1.20431738098716E+04 ,
  84.                   -2.06482942053253E+04 ,
  85.                   -1.50863022876672E+04 ,
  86.                   -1.51383183411507E+03 ,
  87.  
  88.                   -1.03770165173298E+04 ,
  89.                   -9.82710228142049E+05 ,
  90.                   -1.97183011586092E+07 ,
  91.                   -8.73167543823839E+07 ,
  92.                    1.11938535429986E+08 ,
  93.                    4.81807710277363E+08 ,
  94.                   -2.44832176903288E+08 ,
  95.                   -2.40798698017337E+08 ,
  96.  
  97.                    8.06588089900001E-04 ,
  98.                   -5.94997310888900E-04 ,
  99.                    7.93650067542790E-04 ,
  100.                   -2.77777777688189E-03 ,
  101.                    8.33333333333330E-02   );
  102.  
  103.    Q           : ARRAY [ 1 .. 24 ] OF REAL =
  104.                  ( 1.00000000000000E+00 ,
  105.                    4.56467718758591E+01 ,
  106.                    3.77837248482394E+02 ,
  107.                    9.51323597679706E+02 ,
  108.                    8.46075536202078E+02 ,
  109.                    2.62308347026946E+02 ,
  110.                    2.44351966250631E+01 ,
  111.                    4.09779292109262E-01 ,
  112.  
  113.                    1.00000000000000E+00 ,
  114.                    1.28909318901296E+02 ,
  115.                    3.03990304143943E+03 ,
  116.                    2.20295621441566E+04 ,
  117.                    5.71202553960250E+04 ,
  118.                    5.26228638384119E+04 ,
  119.                    1.44020903717009E+04 ,
  120.                    6.98327414057351E+02 ,
  121.  
  122.                    1.00000000000000E+00 ,
  123.                   -2.01527519550048E+03 ,
  124.                   -3.11406284734067E+05 ,
  125.                   -1.04857758304994E+07 ,
  126.                   -1.11925411626332E+08 ,
  127.                   -4.04435928291436E+08 ,
  128.                   -4.35370714804374E+08 ,
  129.                   -7.90261111418763E+07   );
  130.  
  131. BEGIN (* ALGama *)
  132.  
  133.  
  134.                                    (* Initialize *)
  135.  
  136.       Algval := Rinf;
  137.       Scale  := 1.0;
  138.       Alinc  := 0.0;
  139.       Frac   := 0.0;
  140.       Rarg   := Arg;
  141.       Iof    := 1;
  142.       Qminus := FALSE;
  143.       Qdoit  := TRUE;
  144.                                    (* Adjust for negative argument *)
  145.  
  146.       IF( Rarg < 0.0 ) THEN
  147.          BEGIN
  148.  
  149.             Qminus := TRUE;
  150.             Rarg   := -Rarg;
  151.             Top    := Int( Rarg );
  152.             Bot    := 1.0;
  153.  
  154.             IF( ( INT( Top / 2.0 ) * 2.0 ) = 0.0 ) THEN Bot := -1.0;
  155.  
  156.             Top    := Rarg - Top;
  157.  
  158.             IF( Top = 0.0 ) THEN
  159.                Qdoit := FALSE
  160.             ELSE
  161.                BEGIN
  162.                   Frac   := Bot * PI / SIN( Top * PI );
  163.                   Rarg   := Rarg + 1.0;
  164.                   Frac   := LN( ABS( Frac ) );
  165.                END;
  166.  
  167.          END;
  168.                                    (* Choose approximation interval *)
  169.                                    (* based upon argument range     *)
  170.  
  171.       IF( Rarg = 0.0 ) THEN
  172.          Qdoit := FALSE
  173.  
  174.       ELSE IF( Rarg <= 0.5 ) THEN
  175.          BEGIN
  176.  
  177.             Alinc  := -LN( Rarg );
  178.             Scale  := Rarg;
  179.             Rarg   := Rarg + 1.0;
  180.  
  181.             IF( Scale < Zeta ) THEN
  182.                BEGIN
  183.                   Algval := Alinc;
  184.                   Qdoit  := FALSE;
  185.                END;
  186.  
  187.          END
  188.  
  189.       ELSE IF ( Rarg <= 1.5 ) THEN
  190.          Scale := Rarg - 1.0
  191.  
  192.       ELSE IF( Rarg <= 4.0 ) THEN
  193.          BEGIN
  194.             Scale  := Rarg - 2.0;
  195.             Iof    := 9;
  196.          END
  197.  
  198.       ELSE IF( Rarg <= 12.0 ) THEN
  199.          Iof := 17
  200.  
  201.       ELSE IF( Rarg <= RMAX   ) THEN
  202.          BEGIN
  203.  
  204.             Alinc  := ( Rarg - 0.5 ) * LN( Rarg ) - Rarg + Xln2sp;
  205.             Scale  := 1.0 / Rarg;
  206.             Rarg   := Scale * Scale;
  207.  
  208.             Top    := P[ 25 ];
  209.  
  210.             FOR I := 26 TO 29 DO
  211.                Top    := Top * Rarg + P[ I ];
  212.  
  213.             Algval := Scale * Top + Alinc;
  214.  
  215.             Qdoit  := FALSE;
  216.  
  217.          END;
  218.  
  219.                                    (* Common evaluation code for Arg <= 12. *)
  220.                                    (* Horner's method is used, which seems  *)
  221.                                    (* to give better accuracy than          *)
  222.                                    (* continued fractions.                  *)
  223.  
  224.       IF Qdoit THEN
  225.          BEGIN
  226.  
  227.             Ilo    := Iof + 1;
  228.             Ihi    := Iof + 7;
  229.  
  230.             Top    := P[ Iof ];
  231.             Bot    := Q[ Iof ];
  232.  
  233.             FOR I := Ilo TO Ihi DO
  234.                BEGIN
  235.                   Top    := Top * Rarg + P[ I ];
  236.                   Bot    := Bot * Rarg + Q[ I ];
  237.                END;
  238.  
  239.             Algval := Scale * ( Top / Bot ) + Alinc;
  240.  
  241.          END;
  242.  
  243.       IF( Qminus ) THEN Algval := Frac - Algval;
  244.  
  245.       ALGama := Algval;
  246.  
  247. END   (* ALGama *);
  248.